#include <vector>
#include <string>
#include <set>
#include <iostream>
#include <cstdio>
#include <cassert>
#include <numeric>
#include <algorithm>
#include <complex>
#include <ctime>
#include <fstream>

#include <src/indexing.h>
#include <src/configuration.h>
#include <src/lanczos.h>
#include <src/hamiltonian_actor.h>
#include <src/numbers.h>
#include <src/random.h>

int main(int argc, char ** argv) {

  // set random seed
  pcps::set_seed( std::time(NULL) );

  // check to make sure the correct number of arguments were given
  if (argc != 2)
    throw pcps::Exception("Wrong number of arguments.  Correct usage:  eds.exe input_file");

  // create some variables that will be set by the input file
  int nthread, n, cps, t, max_iter, restart_freq, nbond, ntri, n_tp_vecs, tp_tri_num;
  double total_sz;
  bool restart;
  std::vector<int> bonds;
  std::vector<double> bond_strengths;
  std::vector<int> triangles;

  // read input file
  {

    // open the input file
    std::ifstream input(argv[1]);
    if (!input.is_open())
      throw pcps::Exception("Failed to open input file \"%s\".") % argv[1];

    // create a string for reading unused parts of the input file
    std::string not_used;

    // read number of threads
    input >> not_used >> nthread;
    if ( nthread < 0 || nthread > 8 )
      throw pcps::Exception("nthread must be in the range [0, 1, 2, ..., 8]");

    // read whether to restart from file on first lanczos run
    input >> not_used >> restart;

    // read number of sites
    input >> not_used >> n;

    // read configurations per site
    input >> not_used >> cps;

    // read total sz
    input >> not_used >> total_sz;

    // compute total occupation
    t = int( total_sz + 0.5 * n * ( cps - 1 ) + 0.5 );

    // read maximum iteration number
    input >> not_used >> max_iter;

    // read restart_frequency
    input >> not_used >> restart_freq;

    // read n_tp_vecs
    input >> not_used >> n_tp_vecs;

    // read tp_tri_num
    input >> not_used >> tp_tri_num;

    // read number of bonds
    input >> not_used >> nbond;

    // resize bond vectors
    bonds.resize( 2 * nbond );
    bond_strengths.resize( nbond );

    // read bonds
    for (int i = 0; i < nbond; i++)
      input >> bonds.at(2*i) >> bonds.at(2*i+1) >> bond_strengths.at(i);

    // read number of triangles
    input >> not_used >> ntri;

    // resize triangle vectors
    triangles.resize( 3 * ntri );

    // read triangles
    for (int i = 0; i < ntri; i++)
      input >> triangles.at(3*i) >> triangles.at(3*i+1) >> triangles.at(3*i+2);

    // close the file
    input.close();

  }

  // print run parameters
  std::cout << std::endl;
  std::cout << boost::format("nthread        %6i")   % nthread       << std::endl;
  std::cout << boost::format("restart        %6i")   % int(restart)  << std::endl;
  std::cout << boost::format("n              %6i")   % n             << std::endl;
  std::cout << boost::format("cps            %6i")   % cps           << std::endl;
  std::cout << boost::format("total_sz       %6.1f") % total_sz      << std::endl;
  std::cout << boost::format("max_iter       %6i")   % max_iter      << std::endl;
  std::cout << boost::format("restart_freq   %6i")   % restart_freq  << std::endl;
  std::cout << boost::format("n_tp_vecs      %6i")   % n_tp_vecs     << std::endl;
  std::cout << boost::format("tp_tri_num     %6i")   % tp_tri_num    << std::endl;
  std::cout << boost::format("nbond          %6i")   % nbond         << std::endl;
  for (int i = 0; i < nbond; i++)
    std::cout << boost::format("                     %6i  %6i  %20.12f")
                 % bonds.at(2*i)
                 % bonds.at(2*i+1)
                 % bond_strengths.at(i)
              << std::endl;
  std::cout << boost::format("ntri       %6i")   % ntri     << std::endl;
  for (int i = 0; i < ntri; i++)
    std::cout << boost::format("                     %6i  %6i  %6i")
                 % triangles.at(3*i)
                 % triangles.at(3*i+1)
                 % triangles.at(3*i+2)
              << std::endl;
  std::cout << std::endl;

  //const int   n = 12;
  //const int cps =  5;
  //const int   t = 24; //(n*(cps-1))/2;

  // get some integer work space
  std::vector<int>     work(  5*cps + 2             , 0 );

  // run lanczos
  if ( n_tp_vecs == 1 ) {

    double lambda;

    const int vec_len = eds::n_integer_solutions(n, t, cps-1, &work.at(0));

    std::vector<std::complex<double> >      x(vec_len, pcps::zero<std::complex<double> >());
    std::vector<std::complex<double> >     v0(vec_len, pcps::zero<std::complex<double> >());
    std::vector<std::complex<double> >     v1(vec_len, pcps::zero<std::complex<double> >());
    std::vector<std::complex<double> >  t_vec(max_iter*max_iter, pcps::zero<std::complex<double> >());
    std::vector<std::complex<double> >  alpha(max_iter, pcps::zero<std::complex<double> >());
    std::vector<std::complex<double> >   beta(max_iter, pcps::zero<std::complex<double> >());
    std::vector<double>                     w(max_iter, 0.0);
    std::vector<std::complex<double> >  xwork(5*max_iter, pcps::zero<std::complex<double> >());
    std::vector<double>                 rwork(3*max_iter, 0.0);

    eds::HamiltonianActor<std::complex<double> > hamc(n, t, cps, bonds, bond_strengths, triangles, nthread);

    bool converged = false;
    for (int total_iter = 0; total_iter < max_iter && !converged; total_iter += restart_freq) {

      const int loop_iter = ( total_iter + restart_freq >= max_iter ? max_iter - total_iter : restart_freq );

      converged = pcps::lanczos(vec_len, loop_iter, lambda, &x.at(0), &v0.at(0), &v1.at(0), &t_vec.at(0),
                                &alpha.at(0), &beta.at(0), &w.at(0), &xwork.at(0), &rwork.at(0), &hamc, restart);

      restart = true;

    }

    hamc.triple_products(vec_len, 1, tp_tri_num);

  // diagonalize triple products
  } else if ( n_tp_vecs > 1 ) {

    const int vec_len = eds::n_integer_solutions(n, t, cps-1, &work.at(0));
    eds::HamiltonianActor<std::complex<double> > hamc(n, t, cps, bonds, bond_strengths, triangles, nthread);
    hamc.triple_products(vec_len, n_tp_vecs, tp_tri_num);

  } else {
    throw pcps::Exception("user entered %i for n_tp_vecs, which should be >= 1") % n_tp_vecs;
  }

//  // test looping over a configuration
//  if (true) {
//
//    std::vector<int>  z_array(  eds::z_array_size(n)  , 0 );
//    std::vector<int> z_lookup(  (n+1)*(n+1)           , 0 );
//    std::vector<int>  offsets(  (t+1)*(n+1)*(n+1)*cps , 0 );
//    std::vector<int>  spacing(  (t+1)*(n+1)*cps       , 0 );
//
//    eds::prepare_z_array(n, &z_array.at(0), &z_lookup.at(0));
//    eds::prepare_spacing_and_offsets(n, cps, t, &offsets.at(0), &spacing.at(0), &work.at(0));
//
//    std::cout << std::endl;
//    const int correct_total = eds::n_integer_solutions(n, t, cps-1, &work.at(0));
//    printf("  %10i     = correct total\n", correct_total);
//    std::cout << std::endl;
//
//    std::vector<int> ind_vec(correct_total, -1);
//
//    int count = 0;
//    for (eds::Configuration config(n, t, cps); !config.finished(); config++) {
//
//      // compute compound index
//      const int cmpd = eds::index_from_occ_vec(n, cps, t, &offsets.at(0), &spacing.at(0),
//                                               &z_lookup.at(0), &z_array.at(0), config.c_ptr(), &work.at(0));
//
//      // store the index
//      ind_vec.at(count) = cmpd;
//
//      //// print compound index
//      //printf("  %10i    ", cmpd);
//
//      //// print configuration
//      //for (int i = 0; i < n; i++)
//      //  printf(" %2i", config[i]);
//      //std::cout << std::endl;
//
//      if (++count % 1000000 == 0) {
//        printf("  %10i", count);
//        std::cout << std::endl;
//      }
//
//    }
//
//    // sort the indices
//    std::sort(ind_vec.begin(), ind_vec.end());
//
//    // check that the indices are correct
//    int order_error = 0;
//    for (int i = 1; i < ind_vec.size(); i++) {
//      const int x = ind_vec[i] - ind_vec[i-1] - 1;
//      order_error += x*x;
//    }
//
//    std::cout << std::endl;
//    printf("  %10i     = correct total\n", correct_total);
//    printf("  %10i     = actual total\n", count);
//    printf("  %10i     = highest+1\n", *ind_vec.rbegin()+1);
//    printf("  %10i     = lowest\n", *ind_vec.begin());
//    printf("  %10i     = order error\n", order_error);
//    std::cout << std::endl;
//
//    assert( count == correct_total );
//    assert( *ind_vec.rbegin() + 1 == correct_total );
//    assert( *ind_vec.begin() == 0 );
//
//  }
//
//  // test indexing of occupation number vector
//  if (false) {
//
//    std::vector<int> occ_vec(n, 0);
//
//    std::set<int> ind_set;
//
//    for (int i0 = 0; i0 < cps; i0++)
//    for (int i1 = 0; i1 < cps; i1++)
//    for (int i2 = 0; i2 < cps; i2++)
//    for (int i3 = 0; i3 < cps; i3++)
//    for (int i4 = 0; i4 < cps; i4++)
//    for (int i5 = 0; i5 < cps; i5++)
//    //for (int i6 = 0; i6 < cps; i6++)
//    //for (int i7 = 0; i7 < cps; i7++)
//    {
//
//      occ_vec[0] = i0;
//      occ_vec[1] = i1;
//      occ_vec[2] = i2;
//      occ_vec[3] = i3;
//      occ_vec[4] = i4;
//      occ_vec[5] = i5;
//      //occ_vec[6] = i6;
//      //occ_vec[7] = i7;
//
//      if ( std::accumulate(occ_vec.begin(), occ_vec.end(), 0) != t ) continue;
//
//      const int cmpd = eds::index_from_occ_vec(n, cps, t, &offsets.at(0), &spacing.at(0),
//                                               &z_lookup.at(0), &z_array.at(0), &occ_vec.at(0), &work.at(0));
//
//      // print config and index
//      printf("  %10i    ", cmpd);
//      for (int i = 0; i < n; i++)
//        std::cout << occ_vec[i];
//      std::cout << std::endl;
//      //std::cout << std::endl;
//
//      assert( ind_set.count(cmpd) == 0 );
//      ind_set.insert(cmpd);
//
//    }
//
//    assert( ind_set.size() == eds::n_integer_solutions(n, t, cps-1, &work.at(0)) );
//    assert( *ind_set.begin() == 0 );
//    assert( *ind_set.rbegin() == eds::n_integer_solutions(n, t, cps-1, &work.at(0)) - 1 );
//
//  }
//
//  // test z array
//  if (false) {
//
//    int bins  = 12;
//    int balls = 12;
//    const int * const z_ptr = &z_array.at(0) + z_lookup[bins*(n+1)+balls];
//    std::vector<int> occ_vec(bins, 0);
//    for (int i = 0; i < balls; i++)
//      occ_vec[i] = 1;
//
//    const int max_count = eds::binom_coeff(bins, balls);
//    int prev_ind = -1;
//    for (int count = 0; count < max_count; count++) {
//
//      // compute index
//      int cmpd = 0;
//      {
//        int k = 0;
//        for (int l = 0; l < bins; l++)
//          if (occ_vec[l] == 1)
//            cmpd += z_ptr[ (k++)*bins + l ];
//      }
//      assert( cmpd == prev_ind + 1 );
//      prev_ind = cmpd;
//
//      // print config and index
//      for (int i = 0; i < bins; i++)
//        std::cout << occ_vec[i];
//      printf("  %10i", cmpd);
//      std::cout << std::endl;
//
//      // move balls
//      bool flag = false;
//      int accumulation = 0;
//      for (int i = bins-1; i >= 0; i--) {
//
//        if (occ_vec[i] == 1)
//          accumulation++;
//
//        if (occ_vec[i] == 0)
//          flag = true;
//
//        if (occ_vec[i] == 1 && flag) {
//          occ_vec[i] = 0;
//          assert(occ_vec[i+1] == 0);
//          for (int j = 1; j <= accumulation; j++)
//            occ_vec[i+j] = 1;
//          i += accumulation + 1;
//          for ( ; i < bins; i++)
//            occ_vec[i] = 0;
//          break;
//        }
//        
//      }
//
//    }
//
//  }

  return 0;

}
